function b=source(term,n,varargin)

% these are for testing use
% f{1}=@(x) -x.*0;
% f{2}=@(x) -x.*0;
% f{3}=@(x) -x.*0;
% f{4}=@(x) -x.*0;
% linearpartofbjk=ones(term);
% constpartofbjk=ones(term);
% constpartofhjk=ones(term);

if(exist('data.mat','file'))
    load data;
else
    disp('Data file reading error!')
end

unknownsposition=(1:n-1)/n;
h=1/n;
coeff1=1/h^2; coeff2=1/2/h;
tmp=zeros(n-1,term);
leftboundary=zeros(1,term);
rightboundary=zeros(1,term);
for ii=1:floor((nargin-1)/2)
    leftboundary(ii)=varargin{2*ii-1};
    rightboundary(ii)=varargin{2*ii};
end

for row=1:term
    tmp(:,row)=f{row}(unknownsposition);
    for column=1:term
        bjk=@(x) linearpartofbjk(row,column)*x+constpartofbjk(row,column);
        hjk=@(x) constpartofhjk(row,column)+x-x;
        tmpbjk=bjk(unknownsposition); %length is n-1
        tmphjk=hjk(unknownsposition);
        tmp(1,row)=tmp(1,row)-coeff1*tmpbjk(1)*leftboundary(column)+coeff2*tmphjk(1)*leftboundary(column);
        tmp(n-1,row)=tmp(n-1,row)-coeff1*tmpbjk(n-1)*rightboundary(column)-coeff2*tmphjk(n-1)*rightboundary(column);
    end
end

b=reshape(tmp,[(n-1)*row,1]);

